Pronalaženje šablona¶
Još jedan od čestih problema sa kojim se susreće bioinformatika je pronalaženje šablona unutar teksta. Pored egzaktnog pronalaženja šablona u jednoj sekvenci, što je slučaj sa KMP algoritmom, ili pronalaženja približnih poravnanja, često je potrebno efikasno istovremeno pronaći (egzaktno ili približno) pojavljivanje podniski u većem broju niski ili pojavljivanje većeg broja niski u jednoj niski. Na primer, realan je problem pronalaženja sekvenci koje kodiraju jedan gen u više uzoraka neke vrste.
Jedan pristup rešavanju ovog problema je konstrukcija prefiksnog stabla (eng. Trie)
Trie¶
Struktura Trie konstruiše prefiksno stablo od niski koje se pretražuju (šablona) a zatim "provlači" nisku kroz stablo i pronalazi u niski sva pojavljivanja šablona od kojih je stablo konstruisano.
Inicijalno, stablo je prazno sa korenom koji predstavlja prazni karakter (''), dok će kasnije čvorovi stabla predstavljati karaktere šablona. Šabloni se jedan po jedan dodaju u stablo. Krećući se od korena, karakter po karakter šablona se dodaje u stablo kao dete tekućeg čvora a zatim se dodavanje nastavlja od novododatog čvora. Ukoliko tekući čvor već ima dete-čvor koji odgovara tekućem karakteru, novi čvor se ne dodaje već se pretraga nastavlja od postojećeg.
Pretraga šablona počinje od korena stabla i prvog karaktera niske, svakim sledećim karakterom niske prati se put od tekućeg čvora krećući se po čvorovima koji odgovaraju tekućim karakterima. Ukoliko iz tekućeg čvora ne postoji dete-čvor sa potrebnim karakterom, pretraga se vraća na početak i postupak se ponavlja od sledećeg karaktera. Ako se tokom pretrage dođe do lista, znači da je pronađen šablon kome odgovara list, pretraga se vraća na početak i počinje pretraga od sledećeg karaktera (u odnosu na početni u prethodnoj pretrazi od korena).
Ova struktura ima ograničenje da ni jedan šablon ne sme biti prefiks nekog drugog šablona.
Primer:
Šabloni: ATC, ACT, TG
Stablo: ''
Dodavanje prvog šablona u stablo:
Šablon ATC
Korak 1: Korak 2: Korak 3:
'' '' ''
| | |
A A A
| |
T T
|
CDodavanje drugog šablona u tekuće stablo:
Šablon: ACT
Korak 4: Korak 5: Korak 6:
'' '' ''
| | |
A A A
| | \ | \
T T C T C
| | | |
C C C TDodavanje trećeg šablona u tekuće stablo:
Šablon: TG
Korak 6: Korak 7:
''-- ''--
| \ | \
A T A T
| \ | \ |
T C T C G
| | | |
C T C TPronalaženje šablona u niski AATCTG:
('[]AATCTG'):
*''--
| \
A T
| \ |
T C G
| |
C T
([A]ATCTG):
''--
| \
*A T
| \ |
T C G
| |
C T
([A]ATCTG) Iz tekućeg čvora ne postoji dete-čvor A, povratak na početak:
*''--
| \
A T
| \ |
T C G
| |
C T
(A[A]TCTG):
''--
| \
*A T
| \ |
T C G
| |
C T
(AA[T]CTG):
''--
| \
A T
| \ |
*T C G
| |
C T
(AAT[C]TG):
''--
| \
A T
| \ |
T C G
| |
*C T
(AAT[C]TG) Pronađen šablon ATC na poziciji 4 - duzina(ATC) = 4 - 3 = 1, povratak na početak:
*''--
| \
A T
| \ |
T C G
| |
C T
(AA[T]CTG):
''--
| \
A *T
| \ |
T C G
| |
C T
(AA[T]CTG) Iz tekućeg čvora ne postoji dete-čvor C, povratak na početak:
*''--
| \
A T
| \ |
T C G
| |
C T
(AAT[C]TG) Iz tekućeg čvora ne postoji dete-čvor C, povratak na početak:
*''--
| \
A T
| \ |
T C G
| |
C T
(AAT[C]TG) Iz tekućeg čvora ne postoji dete-čvor C, povratak na početak:
*''--
| \
A T
| \ |
T C G
| |
C T
(AATC[T]G):
''--
| \
A *T
| \ |
T C G
| |
C T
(AATCT[G]):
''--
| \
A T
| \ |
T C *G
| |
C T
Pronađen šablon TG na poziciji 6 - duzina(TG) = 6 - 2 = 4, nema više karaktera, kraj:Struktura čvora Trie stabla:
import random
class TrieNode:
def __init__(self, label):
self.label = label
self.neighbors = {}
self.is_leaf = True
self.id = random.randrange(10000000) # Opcioni, id čvora
# Dodavanje čvora suseda (dete-čvora) koji predstavlja zadati karakter
def add_neighbor(self, character):
self.is_leaf = False
self.neighbors[character] = TrieNode(self.label + character)
# Provera da li čvor ima suseda sa zadatim karakterom
def has_neighbor(self, character):
if character in self.neighbors:
return True
return False
# Funkcija vraća susedni čvor tekućeg čvora sa zadatim karakterom
def get_neighbor(self, character):
if self.has_neighbor(character):
return self.neighbors[character]
return None
# Funkcija vraća sve susedne čvorove tekućeg čvora
def get_neighbors(self):
return self.neighbors
Struktura Trie koja se sastoji od povezanih čvorova:
import networkx as nx
import matplotlib.pyplot as plt
class Trie:
def __init__(self, patterns):
self.root = self.trie_construction(patterns)
# Konstrukcija stabla od zadatih šablona
def trie_construction(self, patterns):
root = TrieNode('')
for pattern in patterns:
current_node = root
for current_symbol in pattern:
if not current_node.has_neighbor(current_symbol):
current_node.add_neighbor(current_symbol)
current_node = current_node.get_neighbor(current_symbol)
return root
# Funkcija vraća koren stabla
def get_root(self):
return self.root
# Funkcija pretražuje šablone od tekućeg prefiksa (tekuće pozicije u niski)
# i vraća poziciju prvog pronađenog šablona
def prefix_trie_matching(self, text):
current_node = self.root
for c in text:
if current_node.is_leaf:
return current_node.label
if current_node.has_neighbor(c):
current_node = current_node.get_neighbor(c)
else:
return False
if current_node.is_leaf:
return current_node.label
else:
return False
# Pronalaženje šablona u zadatoj niski od svake pozicije.
# Funkcija vraće pronađene šablone zajedno sa pozicijama na kojima su pronađeni
def trie_matching(self, text):
results = []
n = len(text)
for i in range(n):
match = self.prefix_trie_matching(text[i:])
if match:
results.append((match, i))
return results
# (OPCIONO)
# Funkcija iscrtava stablo
def draw(self):
edges = []
visited = set([])
stack = [self.get_root()]
Gx = nx.Graph()
node_labels = {}
while len(stack) > 0:
v = stack[-1]
if len(v.label) > 0:
label = v.label[-1]
else:
label = v.label
node_labels[v.id] = label
neighbors = v.get_neighbors()
for character, neighbor in neighbors.items():
Gx.add_edge(v.id, neighbor.id)
found = False
for character, neighbor in neighbors.items():
if neighbor not in visited:
visited.add(neighbor)
stack.append(neighbor)
found = True
break
if not found:
stack.pop()
nx.draw(Gx, with_labels=True, labels=node_labels, font_weight='bold')
plt.show()
# Primer
text = 'bananas'
patterns = ['ab', 'bananas', 'ananas', 'cucumber']
trie = Trie(patterns)
trie.draw()
result = trie.trie_matching(text)
print(f'Text: {text}')
print(f'Found patterns: {result}')
Još jedna korisna struktura je sufiksni niz (bliski rođak sufiksnog stabla). Ovog puta, umesto konstrukcije strukture od šablona koji se pretražuju, struktura sufiksnog niza biće konstruisana od niske u kojoj se traže šabloni. Sufiksni niz predstavlja niz leksikografski sortiranih sufiksa zadate niske. Kako bi se naznačio kraj niske, kao poslednji karakter se dodaje specijalni karakter koji se ne pojavljuje u regularnoj niski. U narednim primerima biće korišćen karakter $.
Pronalaženje šablona pomoću sufiksnog niza vrši se binarnom pretragom, tražeži sufikse čiji je prefiks jednak traženom šablonu. Pri pronalaženju jednog traženog sufiksa, pretraga se, linearno, proširuje na susedne sufikse u nizu čiji prefiksi takođe mogu odgovarati šablonu koji se pretražuje.
Moguće je konstruisati sufiksni niz pomoću više niski, pri čemu je moguća istovremena pretraga šablona u više niski.
class Suffix_array:
def __init__(self, text):
text = text + '$'
self.suffix_array = self.generate_suffix_array(text)
# Izdvajanje svih sufiksa zadate niske
def generate_suffixes(self, text):
suffixes = []
n = len(text)
for i in range(n):
suffixes.append(text[i:])
return suffixes
# Generisanje sufiksnog niza leksikografskim sortiranjem
# niza svih sufiksa zadate niske
def generate_suffix_array(self, text):
raw_suffixes = self.generate_suffixes(text)
n = len(text)
suffixes = list(zip(raw_suffixes, [i for i in range(n)]))
suffixes.sort()
return suffixes
# Funkcija proverava da li je niska 1 prefiks niske 2
def is_prefix(self, text_1, text_2):
n = len(text_1)
m = len(text_2)
if n > m:
return False
return text_1 == text_2[:n]
# Pronalaženje pojedinačnog šablona pomoću sufiksnog niza
def pattern_matching_with_suffix_array(self, pattern):
n = len(self.suffix_array)
l = 0
d = n
while l <= d:
mid = (l + d) // 2
mid_suffix = self.suffix_array[mid][0]
if self.is_prefix(pattern, mid_suffix):
i = mid - 1
j = mid + 1
while i >= 0:
if not self.is_prefix(pattern, self.suffix_array[i][0]):
i += 1
break
i -= 1
while j < n:
if not self.is_prefix(pattern, self.suffix_array[j][0]):
break
j += 1
return [x[1] for x in self.suffix_array[i:j]]
if mid_suffix < pattern:
l = mid
else:
d = mid
return []
# Primer
text = 'panamabananas'
pattern = 'an'
arr = Suffix_array(text)
positions = arr.pattern_matching_with_suffix_array(pattern)
print(f'Text: {text}')
print(f'Suffix array: ')
for suffix, position in arr.suffix_array:
print(f'\t{suffix}: {position}')
print(f'Pattern: {pattern}')
print(f'Positions of pattern in text: {positions}')
Kako nije najefikasnije čuvati ceo sufiksni niz (kvadratna prostorna složenost po broju karaktera u niski), pametnom reprezentacijom se niz može sačuvati u efikasnijem obliku i rekonstruisati tokom pretrage.
Transformacija koja omogućava ovakvu uštedu je Burrows-Wheeler transformacija.
Početni korak transformacije je generisanje svih cikličnih permutacija ulazne niske
Primer:
Niska: ABC$
Permutacije:
ABC$
$ABC
C$AB
BC$ALeksikografskim sortiranjem dobijenih niski dobija se matrica
$ABC
ABC$
BC$A
C$ABOva matrica odgovara sufiksnom nizu:
$
ABC$
BC$
C$Poslednja kolona dobijene matrice dovoljna je za rekonstrukciju polazne niske i sufiksnog niza i predstavlja rezultat BW transformacije:
C
$
A
BRekonstrukcija sufiksnog niza:
BW:
C
$
A
BSortiranjem se dobija prva kolona matrice
$ C
A $
B A
C BPoslednji karakter niske je $. U poslednjoj koloni $ se nalazi u drugom redu, kako se karakter A nalazi u prvoj koloni istog reda, sećajući se postupka formiranja table od cikličnih permutacije niske, taj karakter A sledi nakon karaktera $ u cikličnoj permutaciji. Zato je:
$A C
A $
B A
C BKarakter A se u poslednjoj koloni nalazi u trećem redu. Prvi karakter u istom redu odgovara karakteru koji sledi nakon karaktera A u cikličnoj permutaciji. To je, u ovom slučaju, karakter B:
$ABC
A $
B A
C BTime je polazna niska rekonstruisana
BWT se može koristi za pretragu teksta, baš kao što je slučaj sa sufiksnim nizom. U prethodnom primeru se karakter A nalazi u tekstu samo jednom. Ipak, može se dogoditi da se jedan karakter pojavljuje više puta u niski (što najčešće i jeste slučaj) pa je potrebno znati kom po redu karakteru odgovara karakter iz poslednje kolone.
Primer:
Tekst: ABA$
ABA$
$ABA
A$AB
BA$A
$ A
A B
A $
B APrvom karakteru A iz poslednje kolone odgovara prvi karakter A iz prve kolone. U opštem slučaju, $i$-ti karakter $c$ iz poslednje kolone odgovara $i$-tom istom karakteru iz prve kolone.
Ovo svojstvo mapiranja poslednje kolone na prvu korisno je za pretragu šablona koja se izvodi na način analogan sufiksnom nizu. U nastavku će biti predstavljen algoritam koji broji pojavljivanja traženog uzorka u tekstu.
class BWT:
def __init__(self, text):
text = text + '$'
self.bwt_text = self.bw_construction(text)
first_column = (list(self.bwt_text[:]))
first_column.sort()
self.first_column = ''.join(first_column)
# Generisanje cikličnih permutacija zadate niske
def generate_cyclic_permutations(self, text):
permutations = []
n = len(text)
for i in range(n):
prefix = text[:(n-i)]
suffix = text[(n-i):]
permutations.append(suffix + prefix)
return permutations
# Konstrukcija BW transformisanog teksta
def bw_construction(self, text):
permutations = self.generate_cyclic_permutations(text)
permutations.sort()
bwt_text = ''.join([x[-1] for x in permutations])
return bwt_text
# Inverzna BW transformacija
def inverse_bwt(self):
bwt_text = self.bwt_text
n = len(bwt_text)
last_column = list(bwt_text)
columns = last_column[:]
columns.sort()
result_index = bwt_text.index('$')
for i in range(n - 1):
for j in range(n):
columns[j] = last_column[j] + columns[j]
columns.sort()
return columns[result_index]
# Mapiranje i-tog (index) karaktera poslednje kolone na odgovarajući karakter
# iz prve kolone matrice BW transformacije
def last_to_first(self,index):
first_column = self.first_column
last_column = self.bwt_text
char_at_index = last_column[index]
n = len(last_column)
rank = 0;
for i in range(index+1):
if last_column[i] == char_at_index:
rank += 1
first_column_rank = 0
for i in range(n):
if first_column[i] == char_at_index:
first_column_rank += 1
if first_column_rank == rank:
return i
# Pretraga šablona pomoću BW transformacije
def bw_matching(self, pattern):
first_column = self.first_column
last_column = self.bwt_text
top = 0
bottom = len(first_column) - 1 # Broj redova matrice
m = len(pattern)
j = 0
# Binarna pretraga (L = top, D = bottom)
while top < bottom:
if j < m:
symbol = pattern[j]
j += 1
# Oznaka koje signaliziraju
# da granice novog intervala još nisu postavljene
top_index = -1
bottom_index = -1
for i in range(top,bottom + 1):
# ako je pronađen j-ti karakter šablona u poslednjih koloni,
# evidentira se početna granica novog intervala
if top_index == -1 and last_column[i] == symbol:
top_index = i
bottom_index = i
# evidentira se završna granica novog intervala
elif last_column[i] == symbol:
bottom_index = i
if top_index == -1 or bottom_index == -1:
return 0
# Mapira se interval na interval u prvoj kolonu
# i nastavlja pretraga
top = self.last_to_first(top_index)
bottom = self.last_to_first(bottom_index)
else: # sve dok nije obrađen poslednji karakter
break
# Širina poslednjeg pronađenog intervala se vraća
# i predstavlja broj pojavljivanja traženog šablona
return bottom - top + 1
# Primer
text = 'panamabananas'
pattern = 'an'
bwt = BWT(text)
i = 3
print(f'Text: {text}')
print(f'First column: {bwt.first_column}')
print(f'BWT text: {bwt.bwt_text}')
print(f'Last to first: #{i} in the last column -> number #{bwt.last_to_first(i)} in the first')
print(f'Inverse BWT text: {bwt.inverse_bwt()}')
print(f'Pattern "{pattern}" count: {bwt.bw_matching(pattern)}')